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ABSTRACT 


A new formulation of electromagnetic wave scattering by convex, 
two dimensional conducting bodies was introduced by Kriegsman et al 
in ref.[l]. This formulation, called the On Surface Radiation 
Condition (OSRC) approach is based upon an expansion of the 
radiation condition applied directly on the surface of the 
scatterer. Past approaches involved applying a radiation condition 
at some distance from the scatterer in order to achieve a nearly 
reflection free truncation of a finite difference time domain 
lattice. Application of a suitable radiation condition directly on 
the surface of the convex conducting scatterer can lead to 
substantial simplification of the frequency domain integral 
equation for the scattered field, which is reduced to just a line 
integral. In this thesis, the application of the second order 
radiation boundary condition of ref.[l] to concave and re-entrant 
two dimensional bodies is considered. The specific shapes 
considered are slit circular cylinders and thick semi-circular 
cylinders. The formulation is done for the more challenging case of 
Transverse Electric (TE) polarization. The bodies considered here 
are two dimensional in nature and are perfectly conducting. Results 
are presented for the surface current distribution and scattering 
width of these cylinders and comparison is made with known results. 
It was found that the OSRC approach may not be applicable to 
concave structures and also for higher values of kb for convex 
bodies. 
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I.INTRODUCTION 


The On-Surface Radiation Condition (OSRC) approach 
presented here is a medium frequency technique for modeling 
electromagnetic scattering from complex bodies. It is based 
upon an expansion of the radiation condition applied directly 
on the surface of the scatterer. Past techniques involved 
applying a radiation condition at some distance from the 
scatterer in order to achieve a nearly reflection free 
truncation of a finite time domain lattice. Application of a 
suitable radiation condition directly on the surface of the 
scatterer can lead to substantial simplification of the 
frequency domain integral equation for the scattered field, 
which is reduced to just a line integral. The OSRC technique 
converts the usual surface integral equation for the 
scattering problem into either an integration of unknown 
quantities or a simple ordinary differential equation for 
convex two-dimensional objects. For the Transverse Electric 
(TE) case, the problem reduces to solving an ordinary 
differential equation around the scatterer contour and for the 
Transverse Magnetic (TM) case, the solution is known 
explicitly. It is currently applicable to convex conducting 
cylinders of arbitrary cross sections for both TE and TM 
cases. The OSRC approach was motivated, as mentioned in 
ref.[l], by numerical experiments conducted over the past 


1 






twenty years aimed at simulating scalar or vector wave 
propagation and scattering using a finite difference time 
domain (FD-TD) model of the governing wave equation. 

In the subsequent sections of this thesis, we will use the 
OSRC theory to study scattering by two dimensional, convex and 
concave shaped conducting scatterers for the TE case. The 
specific shapes considered are the slit circular cylinder and 
the thick semi-circular cylinder. More emphasis is placed on 
the TE case using the second order boundary condition. Results 
and comparisons are shown for the TE case. Comparisons of 
current distributions for circular cylinders with nomalized 
circumferences (kb) of 5 and 10 with ref.[l] are made. Also 
comparisons of bistatic cross sections for circular cylinders 
are made for several values of Kb. For the thick semi-circular 
cylinder and the slit cylinders, comparison is shown for the 
induced current density and the normalized scattering width. 
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II. SCATTERING FROM CIRCULAR CYLINDERS: OSRC APPROACH 


In this chapter we will first apply the OSRC technique to 
two dimensional, perfectly conducting and convex shaped 
cylinders. Next we will apply the second order boundary 
conditions to circular cylinders with radius r^^ and obtain a 
second order differential equation for the surface scattered 
field. The differential equation is solved by the Fourier 
series method. Finally expressions are obtained for the 
surface current distribution and the scattered fields. 

A. GENERAL FORMULATION FOR THE TE POLARIZATION CASE 

This formulation is based on the second order radiation 
conditions mentioned in ref.[l]. For the TE polerization case, 
the problem reduces to solving an ordinary differential 
equation around the scatterer surface contour. 

Consider a plane electromagnetic wave illuminating a two 
dimensional, perfectly conducting, convex shaped cylinder for 
the TE polarization case. The incident wave, propagating at an 
angle a with respect to the -x axis (see Fig. 1), is given by 

ff =[/ • [/ (1) 

"inc ^ ' 
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Fig.l TE incident field 

where z is the unit vector parallel to the cylinder axis; o is 
the frequency of the incident wave (k= u/c); c is the is the 
speed of light in free space; x and y are the corresponding 
cartesian coordinates in the plane orthogonal to 2. 

Since the problem is two dimensional, the scattered 


magnetic field (see Fig.2) is assumed to be 









( 2 ) 


The scattered field at any distance r is known from the 
values of the field U^(T') and its normal derivative 
dUsCr'J/dv * on the surface of the scatterer as (ref. [6]) 

^ dv' 8v' 


y 



Fig.2 Geometry of the problem 

where C represents the boundary of the cylinder's cross 
section; 0/dv' denotes an outward normal derivative on C; G is 
the two dimensional free space Green's function, 
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(4) 


Gir/D 




ikR) 


R=\F-z‘\=^J{x-x') *+ {y-y')^ < 5 ) 

T and r’ are co-ordinates of the field and source points, i.e. 
^=(x,y) and P = (x', Y') » and s* is the arc length parameter. 

For the TE case, vanishing of the total tangential 
electric field on the surface of the cylinder produces 

dv' 0v^ 

The objective is to determine Uj on the surface of the 
scatterer. For later purposes, it is necessary to obtain an 
expression for the far-zone fields. Using the asymptotic form 
of the Hankel function it can be shown that the far-zone 
scattered field is 


with 


Us 



v/lTcW ^ c dv 


(7) 


( 8 ) 


where 

cosh=V.f and cosi|r= (H.. 

The true fields satisfy Sommerfeld’s radiation condition 
(ref.[6]) on an infinitely large outer cylinder. However, in 
the OSRC technique, a higher order boundary condition such as 
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B 2 (defined below) is applied directly on the surface of the 
object. Implicit in this approximation is the fact that only 
outgoing waves are assumed to exist on the surface. Thus, the 
approximation is expected to have serious drawbacks in 
geometries where significant stored energy can exist such as 
in re-entrant and slit cylinders. In the sections to follow, 
we will be dealing precisely with these situations. 

The form of second order boundary condition Bj is given 
from (ref.[l]) with -jk in place of jk as 

C(s') , 1 y_a .g. 

2 B[jk+C(s')] 2[jk+C(sO] dv' 

where C(s') is the curvature of the surface at s'. 

The condition B 2 U 5 = 0 is applied to circularly shaped 

cylinders in the next section. 

B. APPLICATION TO CIRCULAR CYLINDERS 

At the surface of the perfect electric conductor, the 
boundary conditions on the total fields state that 

iixH=Js ; j?xF=0 (10) 

where Jj is the total surface current. For the TE case the 
scattered fields are 
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( 11 ) 


h:=u^ 


and 



1 dUs 

JUHo dv' 


Hence we have 9 (U,^+U5)/3v'=0 and J^=“(U 5 +U,^) on the 
surface, where v' is positive for convex surfaces and negative 
for concave surfaces. 

For a=0 in (1), the normal derivative of Uj can be 
rewritten as 

^ ~-h =±jkcos^'U^„^ . ( 12 ) 

9v' (±)dr' dr' 


Substituting the expression for and (12) in (9) , setting 
the curvature C(s') equal to 1/r^, and replacing by 


gives 


where 




(13) 


C, = ±- 


Sk^rl-Bjkr, and 


Bk^rl-12jkZg-3 


C,=- 


* 8ic*ro-12j)cr„-3 


(14) 


In C,, the upper sign is valid for convex surfaces and the 
lower one for concave surfaces. The above equation is valid on 
the surface r*=r^ with being 2n periodic. Equation (13) is 
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linear with constant coefficients and can be solved using 
standard analytical or numerical methods. Analytical methods 
such as the method of undetermined coefficients and variation 
of parameters yielded complicated integrals for the solution 
of the second order non-homogeneous differential equation 
(13). Therefore (13) was solved using a Fourier series 
expansion. 

1. SOLUTION OF THE DIFFERENTIAL EQUATION FOR THE 
SURFACE SCATTERED FIELD 

Rewriting (11) in terms of t/g gives 

^2-^-t^s=C'iCOs4)^t7,„^ • 

The incident field on the surface of the cylinder is 


^inclr'.r, 


-Q-Jkr'coBif 


(16) 


which may be expressed in a Fourier series as (ref.[6]) 

■ ‘I’* 

We will assume a complex exponential form for the Fourier 
series representation of (7g, 
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(18) 


Us= E 


a— 


The coefficient will be solved by substituting into 
equation (15). For Ug on the circular portion of the cylinder, 
we have 


dU. 

15 


7=jE 


(19) 


and 


a^- ” 


( 20 ) 


Substituting these into the Left Hand Side (LHS) of (15), 
we get 




X2«.« 


( 21 ) 


Now expanding the Right Hand Side (RHS) of (15), 


qc/,„,cos4)'=qf; j-"J-„(/cr,) (e^’+'+e'^*') 




We now let n+l=in -» n=in-l in the first series and n-l=in 
n=in+l in the second to obtain 


10 










(23) 


In obtaining (23) , use was made of the expression 
. Equating the coefficients of (21) and (23) 


yields the following expression for the Fourier series 
coefficient A^, 






l + Cjil* 


(24) 


We note that the coefficients are dependent on the radius r^, 
of the cylindrical surface. For numerical calculations, the 
series in (18) needs to be truncated. For determining the 
truncation value of n, we determine A„ for large n. For large 
values of the order n, we have from ref.flO] 

J„{kr„) - _i:_(.f^)". (25) 

2n 


Also using the recurrence relation 


1 ' 


(26) 


we obtain for 


ekr^ 

~2ir 


< 1 
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1 


(27) 


2^2icn 2/j 


(lh£o)n-l 


Therefore for large n, can be expressed as 




2,/2x 2 


n 


(28) 


Hence decays very rapidly, for example 


!io _ 3^0-11 


2. SURFACE CURRENT DISTRIBUTION 

Having found the surface magnetic field, it is a simple 
task to determine the surface current. The total electric 
surface current for the circular cylinder is given by 

■ (29) 

Substituting (17) and (18) in (29) gives 

«• 

. (30) 

After expanding and simplifying, (30) reduces to 
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{krj 1 cosiKb' 


(31) 



where e„ = 1 for n equal to zero and c„ = 2 for n not equal to 
zero, and C, and Cj are given in (14) . However for the 
positive sign is taken since the surface is convex. The 
magnitude of the normalized surface current can be expressed 
as 


-7^ 1= |-£ i 1 COSn^'\ . 

tilne 0 l + 


3. FAR ZONE SCATTERED FIELDS 

The far-zone scattered magnetic field for the TE 
case given at (8) can also be expressed as (ref.[5] page 374) 

^ ^ 4 nk Jc dv' 

where C is a closed contour around the cylinder. 

In our case, with reference to Fig.2, we have 


13 








9 


(34) 


gilc/.r^_ gi*:^.co8*_gi*r,co« (V-*) 




(35) 


(jp. O') =±COS (<j>'-(J)) , 


(36) 


dlil 

av' 




(37) 


Substituting (34) through (37) in (33) and simplifying we 

In the derivation of the above equation, we made use of 

the fact that yicos((t'-<l»)A*“'*''*>»-i-U, . 

dr' * 


For a full cylinder we can further simplify equation (38) 
by solving the integral portion as 

2 R —JnM 

/e''*^-”’(*'-*’e^'^d(|)'=(-l)"2iie » eJ''*J„{kr„) . (39) 

0 

Substituting the above value in (38) and after some 
manipulation, we get 
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(40) 




+j ^A„y„{kr^)]co3n^ . 


We list below a summary of the various quantities that 
were determined in this chapter. 


Surface scattered field. . .U, 




/=r -■ 


Slr'=r, 




Total surface field. . • 1 

•• 

Surface current . ‘^♦“"52 [Aj^+j'^Jj^ikr^)] e^"*' 


Fourier coefficient . A„= 


CJ-(^*^^ji^(kr,) 




- 


1 ^ i 


1 I \ n-1 / _±_ \ 

2y/2n 2n ' ^n 


s , e'i’^ 


Far fields. . 


/r 




jktoCOB (♦'-♦) 


k:“ dr' Jc 


In the next chapter we will develop expressions for the 
case when the cylinder is comprised of more than one circular 
surface. 
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III. SCATTERING FROM CIRCULAR SLIT CYLINDERS: OSRC APPROACH 


In this chapter we will consider the case of a slit 
cylinder with finite wall thickness. It is assximed that the 
thickness of the wall is electrically small so that the fields 
do not vary along it. Furthermore, we ignore the effect of the 
edges. Fig.3 shows the geometry of the slit circular cylinder. 
On the outer surface C^, the body is convex, while on the 
inner surface Cj it is concave. Since both outward and inward 
travelling waves can exist in the interior region, and since 
the OSRC technique assumes waves of only one kind, it is not 
expected to yield correct results for the geometry. We will 
use the results developed in the previous sections for the 
circular surfaces. To distinguish between the outer and inner 
surface we shall use superscripts (b for the outer radius and 
a for the inner radius) for the Fourier coefficients. 

A. DETERMINATION OF SURFACE CURRENTS 

The expression for the current distribution is given in 
(29) i.e. = -(Uj + Uinc)' where 
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Fig.3 Circular slit cylinder 


In Fig. 3, 0p, and Bpj are the first and second slit angles 
from the -x axis respectively; b and a are the outer and inner 
radii of the cylinder respectively; s is the arc length and is 
taken as zero at <P*=n-B^ at the outer contour C,; the angles 
a and p are considered zero at ^'~»r- 0 p 2 and 
respectively; s=ba and s=a/3 are the arc lengths of the outer 
and inner surfaces respectively. 









For the outer contour C, we get from fig. 3, 




and .e-J"”-'-’--’ 


Similarly for the inner contour we get 


U,.-£ a‘‘’ 


, IJ =«-lkaco«*'| (46) 

and ^inc ® ® • ' ' 

Therefore the current distribution for the outer and inner 
contours can be obtained by substituting (43) and (44) for the 
outer surface, and (45) and (46) for the inner surface 
respectively in (29). After simplification the following 
expression is obtained for the magnitude of the normalized 
surface currents. 

For the outer contour C, we have 

"inc 0 

where a varies from 0 to 2»r-0^,-0^,; c/*” and in are 

given by the expressions from (16) as 


ek^b^-Bjkb ^{b) _4_ j.Q. 

‘ ” Bk^b^-12jkb-3 * " Bk^b^-12jkb-2 
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For the inner contour Cj we have 

^inc 0 

where /? varies from 0 to STr-e^^-e^,; and are similar 

except that a negative sign is taken for C, in (12) since the 
inner contour is a concave surface. 

Along the linear portions Cj and C^, we assume that the 
current varies linearly between the inner and outer curved 
surfaces. 

B. DETERMINATION OF THE FAR FIELDS 

From Fig.4 the total scattered magnetic field for the 
hollow circular slit cylinder with finite thickness will be 


H 






+ 



4 - 



( 50 ) 


Where C, and Cj are the outer and inner contours with radii b 
and a respectively; 0^,-0,j is the total slit extent. 
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Fig.4 Slit cylinder geometry 


From Fig.4 the integration limits in (50) would be 


/c= 4 * 4 * 4 * 4 


(51) 


r^=b-a ro=a 

For the convex contour C,, (38) is used with rp=b and 


A =A where 

n n 


, <W_ 


l + C/*’j3* 


(52) 


Substituting these values in (38) we get 
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( 53 ) 


h! 1 (11) 2j^[J2j-»j^(kb) f <♦'-♦> 

/r 4 nx " Jci 

eJtr'co.<*'-*>. 

X OX' •'Ci 

The integral over is evaluated in the following 

fashion: 


f g JJdbcos (♦'-♦) 

Jci 

2n *“®pl 

_J‘gJlcbcos (♦'-♦) gJn*'- J gJJcbcos (♦'-♦) _ 




( 54 ) 


The first integral can be evaluated in a closed form as 
in the previous chapter. In the second integral we let 
u=7r-0'-*d0'=-du. After some manipulation we get 


r gjkbcoB (♦'-♦) 


Substituting (55) in (53) and carrying out some manipulation 
we get 
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( 56 ) 


lit = -k-^^ 4 (^"i> [f} 2ne^ [ {kb) {kb) + j 


(« 


‘Cl 


/r 4 nk 


.Jn{kb)'\co3n^-Y^e^’*jii{kb) f cos (jiu) du 

+j]^ "f e"-^*^®"‘“**’cos {nu) du] . 


‘’Dl 


For the inner contour Cj with radius a, which has a 
concave surface we have 

-^=k'£j-«4{ka)e^^ , 

N!='f^A„{a)eJ‘^' , 

.(a) Cl^^j-^”*^>y„{ka) 

(/.O') =-cos((j)'-({>) , 
and f ds'=f ad<i>'=-af . 

• Cj ^ Cj V j 

Substituting these values in ( 33 ) gives 


//; =-k ^ 

‘ci 




jkz„COB (V-*) 


Following the same procedure as for the contour C,, we can 
simplify the integral portion and after some manipulations we 
get 
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^ kt * O * ^ ** 

//4=-;c-^ ^ ) 2 a 2nc^ [ (-1) "Jj; ika) J„ (ka) -j 

w ^P* 

.j'„{ka) 1 co3n(|)-52 ”ji(/ca) f e‘-^**®®“*“***cos (nu) du 

e° 

«• ®Pi 

(-1)” fcos (u+(J)) cos (nu) du] . 

° 


( 58 ) 


For solving the linear portions and C^, we can ignore 
the edge effects since the OSRC approach does not currently 
provide edge current singularity behaviour as mentioned in 
ref.[l]. Therefore, for the linear portions we let Ug be equal 
to the average of Uj on the contours C, and Cj at the position 
of the linear portion. 

For contour Cj we have 


Us\c = 


^5lc/^5lc, 




-J 


^ = -J- (e-J*^'‘^«-*') =-jJtcos (7i-e_,) (60) 

dv' dr' 


and (jf. O') =cos (<!)'-({>) =cos (n- 0 p 2 -<|)) . (61) 

Substituting these values in (33) and simplifying we get 
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<=-^4(44) M 


COS0 


P3 


j (coaepj-coB (0,^4^^) ] 


4 nic cos0pa-cos (0p,+4>) 

_gjtMc..e„-c„ («„.♦, 1 ] ] e„ (-1) "cosne 

[g*J*:«co8(epj4*)_^-iWx;oB(ep,4*) j j 


P2 


( 62 ) 


Similarly for the linear portion we have from Fig.4 
U \ +U \ " a <*>4.a <«> 


^ = - jiccos (n -0p,) e(64) 
(/.O')=cos(n-0p,-<}>) . (65) 


Substituting values in (33) and after simplification we get 

J. {A2. \ ■? r__ r i>:«[cos0p,-co8(ep,44.)) 

^ 4 nJc ^ cos0pj-cos (0pj+4>) ^ 

_gJ*McoB0p,-coB(ep,44)Ij [jlL_:^]e„(-l)"cos720p, 

[g-J*«co8(0p,44)_^-J«>co8(0p,44)j j ^ (66) 

In summary, the total scattered magnetic field for the 
hollow circular slit cylinder is 
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_ - ^itr -1 o -f — ** , 

" [-kb['£ 2ne„lJ„(kb)ji,{kb) 

^X* ^ 0 

+ J 1 cosr«J>^- 52e^"ji(Jcb) recos (nu) du 

•• ®p* 

(-l)"f cos (u+<|>) cos (iiu) 

-ica [ j; 2ne„ [ (-1) «cr„ (Ira) (ka) -j(ira) ] cosn«l> 

0 

e„j ''ifnika) f (nu)du 

° 8pi 

n ®Pa 

(-!)"[ cos (u+<t)) e-^***=®*‘“***cos (/ju) du] 

“ Opi 

_ r _ COs9p3 _ Icos6p,-co8 (Bpi+i) J _ iW>Icos0p,-co8 (Bp,*^) J i 

^ cosOp^-cos (0p2+4)) 

- .(*).. (8) 

-E e„ (-1)" I " , ] cosu0p, (b,,.*) _g-iWx:o.(Bp,.*) 3 3 

0 2 

_ r _ COS0pj _ r i*a(coaBp,-co8(Bp,.*))_ i*McoB0p,-coB(Bp,**))3 

^ COS0p_,-COS (0pi+(|)) 

(-i)«[.fiL_J!il_]cosu0p,le-^*«°’‘®*>‘*^’-e'^^^^ (67) 

0 2 


The scattering width a^-o defined as (ref. [6]) 


aa.p=2TCr ■ — 
l^inc|2 


( 68 ) 


For the scattered field given in (49) we have 
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°2-d=-^\-^ tE {kb) Ji, (kb) *j (Ai?) J cosn^' 

8pi 

"E (^) f e-^"«»“<‘'**>cos (nu) du 

0 n 

•• ®pa 

+jE (-1) ” r COS (u+<j)) COS (JIU) 6“^*^®“ *“*♦* du] 

-ka [^2 2Ke„ [ (-l)^JJka)jiika) y„(ka) ] cosucj) 

0 

"E ®/v7 ^^3) f e-^^*®®“<"*^’cos (i3u) du 

0 J 

®P1 

+jE <-l) ”/cos (U+<J)) e^**®®"‘“**>cos (uu) du] 

[ IcoaB^f-coa (Op,*^)} _ ^jkblcoai^f-coa {Bpg*^) ] j 


-[ 


cpsB 


"pi 


£L 


cos0p2-cos (0p2+4>) 

-J^e„ (-!)"[ " ] cosiiBp^ I e-i*«co8 ^-jicbcos (Op,.*) j j 

n 2 

^£1_ r p,i*«lcosep,-co8 ) 1 _ Jict>Icosepi-co8 (Opi**) 1 -i 

. /n . x\ J 


COS0 


^ COS0pj-CO3 (0pi + <l)) 


• (69) 
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IV. NUMERICAL RESULTS AND DISCUSSION 


All analyses were carried out for the TE polarization case 
and when the incident wave is incident from the -x axis, i.e. 
from 0=^180”. In all cases a and b denote the inner and outer 
radii of the cylinders. 

For the circular cylinder. Figs. 5 and 6 show the 
comparisons of the normalized surface current distribution 
between the present formulation (32) and of ref.[l], both 
using the OSRC approach for kb=5 and 10 respectively. The 
results are comparable, checking the validity of our 
formulation. It was shown in ref.[l], that the results 
approach the exact solution. The computer program to compute 
the normalized surface current distribution using the OSRC 
approach is given in Appendix A. Figs. 7,8 and 9 show the 
comparisons of the bistatic widths for the circular cylinder 
between the present formulation and ref.[8] for kb=l,5 and 10 
respectively. Ref.[8] uses the rigorous moment method 
technique for determining the scattering cross sections of 
infinite cylinders of arbitrary geometrical cross section. The 
results again are in good agreement. The computer program to 
compute the scattering width for circular cylinders using the 
OSRC approach is given in Appendix B. Figs.10 and 11 show the 
comparisons of forward (0*0*) and backward (0=180*) scattering 
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widths respectively, versus kb, between the computed and 
ref.[5] (exact solution). It is seen that the results due to 
the OSRC approach show some deviation at higher frequencies 
(i.e. for large kb) in both cases. The computer program to 
compute a/irh versus kb is given in Appendix C. 

For the hollow circular slit cylinder. Fig.12 shows the 
comparison of the normalized surface current distribution for 
a hollow semi-circular cylinder obtained using the present 
formulation (47) and (49) and that in ref.[7]. The problem of 
scattering in ref.[7] is solved by means of an integral 
equation. In Fig.14, the current density, normalized to unity 
for the magnetic field strength of the incident wave, is 
plotted versus perimeter s normalized to unity for the half 
perimeter. The distance variable s is taken equal to zero at 
the point on the center face nearest the source of the 
incident wave. Moving clockwise around the perimeter, s 
increases in value reaching the value s=l at the center face 
furthest from the source. The results again are comparable to 
the exact solution. The computer program to compute (47) and 
(49) for the circular slit cylinder is given at Appendix D. 
Fig.13 shows the comparison of the scattering width. The 
results obtained using the OSRC approach are quite different 
from the exact solutions as shown in Fig.13. This is in spite 
of the good agreement in the magnitude of the current 
distribution. At this point it seemed necessary to carry out 
some phase comparisons for the current distribution, however 
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no known results were available. The computer program in 
Appendix E was thoroughly checked for any bugs and even tested 
for circular cylinders which yielded satisfactory results. It 
is felt that the disagreement in the scattering width due to 
the OSRC method is problably due to the incorrect modeling of 
phase of the current distribution. The OSRC technique assumes 
waves of only one kind (outgoing) to exist on the scattering 
surfaces which may not yield correct results for concave 
structures. 

Fig.14 shows comparisons of forward scattering width 
(0=0") for the hollow circular slit cylinder (b=a and total 
slit angle of 10") between the OSRC approach and the full 
circular cylinder from ref.[5] (exact solution). Fig.15 shows 
the comparison between the OSRC approach and ref.[5] results 
for the hollow circular slit cylinder. Such shapes are 
expected to have cross section resonances, indicating that 
there is interior information contained in the exterior 
scattering data. Here, again the results are quite different. 
Another peculiar result was that we were getting the same 
results for forward and backward scattering widths using the 
computer program given in Appendix F. This might be due to 
some unforseen error in the program. 
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-Cai-unded rcsulls +♦♦ ref.til results 



/?= (degrees) 


Fig.5 Comparison of surface current distribution on a PEC 
circular cylinder with kb=5. 
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Fig. 8 Comparison of bistatic scattering width for the 
circular cylinder with kb=5. 
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Fig. 9 Comparison of bistatla scattering width for the 
circular cylinder with kb=10. 
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-CoMind.ctl results ref.tS] results 



kb 


Fig.10 Comparison of forward scattering width (^=0*) for the 
circular cylinder between the computed (OSRC approach) and 
ref.[5] (exact solution). 
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Fig. 11 Comparison of backward scatting width (^=180*) for the 
circular cylinder between the computed (OSRC approach) and 
ref.[5] (exact solution). 
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Fig.13 comparison of scattering width for the hollow semi¬ 
circular slit cylinder. 









Fig.14 Comparison of forward scattering width for a hollow 
circular slit cylinder with bsa and total slit angle of 10* 
and full circular cylinder from ref.[5] . 






a/nh 


- CoMputed results *** ref#151 results 



Fig.15 Comparison of forward scattering width for the hollow 
circular slit cylinder with b=a and slit angle of 10'. 
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V. CONCLD8ION8 


The OSRC approach was applied to circular and hollow slit 
circular cylinders. For circular cylinders the normalized 
surface current distribution agrees very well with the exact 
solution. Scattering width results, with different values of 
kb were also comparable with the exact solution, as long as kb 
was not large. For the case of hollow circular slit cylinders, 
the OSRC normalized surface current distribution was 
comparable to exact solution results. However, in case of the 
scattering width, the slit cylinder results were quite 
different from the exact solution. This indicated a need to 
carry out phase comparisons of the OSRC surface current 
distribution with the exact solution. 

It is possible that the OSRC approach may not be 
applicable to concave structures and also for higher values of 
kb for convex bodies. As mentioned earlier, both outward and 
inward travelling waves can exist in the interior region of 
the concave scatterer. Since the OSRC approach assumes waves 
of only one kind, it is not expected to be applicable to 
concave bodies. Further the method fails when the scattering 
surface has cavity shaped objects. 






APPENDIX A. COMPUTER PROGRAM TO COMPUTE CURRENT DISTRIBUTION FOR 
CIRCUZAR CYLINDERS 


C THIS PROGRAM COMPUTES THE SURFACE CURRENT DISTRIBUTION FOR A 
C CIRCULAR CYLINDER FOR TRANSVERSE ELECTRIC CASE USING OSRC 
C APPROACH. 

REAL*8 J(0:200),¥(0:200),Dy(0:200),DJ(0:200),K1D,K2D 
REAL*8 Jl(0:200),¥1(0:200),DY1(0:200),DJ1(0:200) 

REAL EN,NMX,P,PHI,K,K1,DTR,PI,FR,B,JPHI 
INTEGER NMAX,N 

COMPLEX CBl,CB2,i,Sl,S2,S3,S4,DN(999),JP 
OPEN(3,FILE='DAT.1•) 

WRITE(*,*) 'ENTER THE WAVE NUMBER K (REAL) IN I/METERS' 
READ(*,*) K 

WRITE(*,*) 'ENTER THE RADIUS (REAL) IN METERS' 

READ(*,*) B 

WRITE(*,*) 'ENTER NMAX (INT)' 

READ(*,*) NMAX 
PI=3.1415927 
DTR=PI/180. 
i=(0,l) 

K1=K*B 

K1D=DBLE(K1) 

CB1=(8.*Kl**2-8.*i*Kl)/(8.*K1**2-12.*i*Kl-3.) 

CB2=-4./(8.*K1**2-12.*i*Kl~3.) 

CALL BES(NMAX,KID,J,Y,DJ,DY) 

C PLOTTING ROUTINE 

WRITE(3,100) 'SURFACE CURRENT DISTRIBUTION FOR CIRCULAR C¥L. ' 

NPHI=91 

PHI1=0. 

PHI2=180. 

WRITE(3,110) NPHI 
WRITE(3,120) PHIl 
WRITE(3,120) PHI2 
DO 33 P=0.,180.,2. 

PHI=DTR*P 

S4=(0,0) 

DO 22 N=1,NMAX+1 
NMX=FLOAT(N)-l. 

IF(NMX.EQ.O.) THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

S1=(0,0) 

S2=(0,0) 

S3=(0,0) 

Sl=CBl*i**(-l) * DJ(NMX)/(1.+CB2*NMX**2) 

S2=J(NMX) 

S3=COS(NMX*(PI+PHI)) 

DN(N)=-(EN*i**(-NMX)*(S1+S2)*S3) 
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S4=S4+DN(N) 

22 CONTINUE 

JP=S4 

JPHI=CABS(JP) 
WRITE(3,120) JPHI 


33 

CONTINUE 


CLOSE(3) 

100 

FORMAT(A) 

110 

FORMAT(15) 

120 

FORMAT(E12.3) 
STOP 


END 
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APPENDIX B. COMPUTER PROGRAM TO COMPUTE a/b FOR CIRCUZAR CYLINDERS 


THIS PROGRAM COMPUTES THE SCATTERING WIDTH FOR A CIRCULAR 
CYLINDER FOR THE TRANSVERSE ELECTRIC POLARIZATION CASE 
USING THE ON - SURFACE RADIATION CONDITION APPROACH. 
ANGLE SI - PI - PHI GOING CLOCKWISE FROM -180 DEGREES 
ALONG THE X-AXIS. 

REAL*8 J(0:200),¥(0:200),DY(0:200),DJ(0:200),K1D 
INTEGER N,NMAX 

REAL TA,TB,K,K1,K2,TP1,TP2,FIL0N,FR,B,A,L,PI,P,NMX,EN,PHI 
COMPLEX DN(999),S2,S3,S4,CN1,i,CBl,CB2 
OPEN(3,FILE*'DAT.2') 

WRITE(*,*) 'ENTER THE WAVE NUMBER(REAL) IN I/METERS' 
READ(*,*) K 

WRITE(*,*) 'ENTER THE OUTER RADIUS(REAL) IN METERS' 
READ(*,*) B 

WRITE(*,*) 'ENTER NMAX(INTEGER)' 

READ(*,*) NMAX 
PI=3.1415927 
DTR=PI/180. 
i=(0,l) 

K1=K*B 
K1D*DBLE(K1) 

CB1=(8.*Kl**2-8.*i*Kl)/(8.*K1**2-12.*i*Kl-3.) 

CB2=-4./(8.*K1**2-12.*i*Kl-3.) 

CALL BES(NMAX,KID,J,Y,DJ,DY) 

C PLOTTING ROUTINE 

WRITE(3,960)'TE CASE PEC CIR. CYLINDER SCATTERING WIDTH ' 

NPHI*91 

PHI1=0. 

PHI2=180. 

WRITE(3,970) NPHI 
WRITE(3,980) PHIl 
WRITE(3,980) PHI2 
DO 50 P =0.,180.,2. 

PHI=DTR*P 

C EVALUATING THE FUNCTION CNl 
S2=(0,0) 

DO 800 N=1,NMAX+1 
NMX=FLOAT(N)-1. 

S3=(0,0) 

S4=(0,0) 

IF(NMX.EQ.O.) THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

S3=J(NMX)*DJ(NMX) 

S4=i**(-1)*CB1*DJ(NMX)**2/(1+CB2*NMX**2) 

DN(N)=2.*PI*EN*(S3+S4)*COS(NMX*(PI-PHI)) 
S2=S2+DN(N) 
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800 CONTINUE 

CN1=S2 

C EVALUATING THE SCATTERING WIDTH DIVIDED BY CYLINDER RADIUS B 
SIGMA=1./(4.*K1)*CABS(-K1*CN1)**2 
WRITE(3,980) SIGMA 


50 

CONTINUE 


CLOSE(3) 

960 

FORMAT(A) 

970 

FORMAT(15) 

980 

FORMAT(E12.3) 


STOP 

END 
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APPENDIX C. COMPUTER PROGRAM TO COMPUTE a/wb versus Icb FOR 
CIRCULAR CYLINDERS 

THIS PROGRAM COMPUTES THE SCATTERING WIDTH FOR A CIRCULAR 
CYLINDER FOR THE TRANSVERSE ELECTRIC POLARIZATION CASE 
USING THE ON - SURFACE RADIATION CONDITION APPROACH. 
(SIGMA /PI*B Vs KB) 

REAL*8 J(0:200),Y(0:200),DY(0:200),DJ(0;200),K1D 
INTEGER N,NMAX 

REAL TA,TB,K,K1,K2,TP1,TP2,FIL0N,FR,B,A,L,PI,P,NMX,EN,PHI,R 
COMPLEX DN(999),S2,S3,S4,CN1,i,CBl,CB2 
OPEN(3,FILE=•DAT.3') 

WRITE(*,*) 'ENTER THE PHI ANGLE'REAL) IN DECREES' 

READ(*,*) P 

WRITE(*,*) ‘ENTER THE OUTER RADIUS(REAL) IN METERS' 
READ(*,*) B 

WRITE(*,*) 'ENTER NMAX(INTEGER)* 

READ(*,*) NMAX 
PI=3.1415927 
DTR=PI/180. 

PHI = DTR*P 
i=(0.1) 

C PLOTTING ROUTINE 

WRITE(3,960) 'TE CASE PEC CIR. CYLINDER SCATTERING WIDTH ' 

NPHI=59 

PHI1=0.1 

PHI2=6. 

WRITE(3,970) NPHI 
WRITE(3,980) PHIl 
WRITE(3,980) PHI2 
DO 50 R =0.1,6.,.1 
K1=R*B 
K1D=DBLE(K1) 

CB1=(8.*Kl**2-8.*i*Kl)/(8.*K1**2-12.*i*Kl-3.) 

CB2=-4./(8.*K1**2-12.*i*Kl-3.) 

CALL BES(NMAX,KID,J,Y,DJ,DY) 

C EVALUATING THE FUNCTION CNl 
S2=(0,0) 

DO 800 N=1,NMAX+1 
NMX=FLOAT(N)-l. 

S3=(0,0) 

S4=(0,0) 

IF(NMX.EQ.0.) THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

S3=J(NMX)*DJ(NMX) 

S4=i**(-1.)*CB1*DJ(NMX)**2/(1.+CB2*NMX**2) 

DN(N)=2.*PI*EN*(S3+S4)*COS(NMX*PHI) 

S2=S2+DN(N) 

800 CONTINUE 
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CN1=S2 

C EVALUATING THE SCATTERING WIDTH DIVIDED BY CYLINDER RADIUS B 
SIGMA«1./(4.*K1*PI)*CABS(-K1*CN1)**2 
WRITE(3,980) SIGMA 
50 CONTINUE 

CLOSE(3) 

960 FOK,'JVT(A) 

970 FORMAT(15) 

980 FORMAT(E12.3) 

STOP 

END 
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APPENDIX D. COMPUTER PROGRAM TO COMPUTE J, FOR SLIT CIRCULAR 
CYLINDERS 

C THIS PROGRAM COMPUTES THE SURFACE CURRENT DISTRIBUTION FOR A 
C CIRCULAR SLIT CYLINDER FOR TRANSVERSE ELECTRIC CASE USING OSRC 
C APPROACH. 

REAL*8 J(0:200),Y(0:200),DJ(0:200),DY(0:200),K1D,K2D 
REAL*8 Jl(0;200),Y1(0:200),DJ1(0:200),DY1(0:200) 

REAL EN,NMX,P,PHI,K,K1,K2,SMAX,PI,FR,A,B,JPHI,L,S,S2,SM2,R 
REAL SMAX2 
INTEGER NMAX,N 

COMPLEX CBl,CB2,CAl,CA2,i,Sl,S3,S4,DN(999),JP,S5 
OPEN(3,FILE='DAT.5') 

WRITE(*,*) 'ENTER THE FREQUENCY (REAL) IN MHZ' 

READ(*,*) FR 

WRITE(*,*) 'ENTER THE FIRST SLIT ANGLE (REAL) IN DEGREES' 
READ(*,*) TA 

WRITE(*,*) 'ENTER THE SECOND SLIT ANGLE (REAL) IN DEGREES' 
READ(*,*) TB 

WRITE(*,*) 'ENTER NMAX (INT)' 

READ(*,*) NMAX 
PI=3.1415927 
i=(0,l) 

DTR=PI/180. 

K=PI*FR/150. 

L=300./FR 
B*.6^L 
A=.5*L 
K1=K*B 
K2=K*A 
K1D=DBLE(K1) 

K2D=DBLE(K2) 

CALL BES(NMAX,KID,J,Y,DJ,DY) 

CALL BES(NMAX,K2D,J1,Y1,DJI,DYl) 

TPl=TA/57.29578 
TP2=TB/57.29578 
CB1=(8.*Kl**2-8.*i*Kl)/(8.*K1**2-12.*i*Kl-3.) 

CB2=-4./(8.*K1**2-12.*i*Kl-3.) 

CA1=-(8.*K2**2-8.*i*K2)/(8.*K2**2-12.*i*K2-3.) 
CA2=-4./(8.*K2**2-12.*i*K2-3.) 

C PLOTTING ROUTINE 

WRITE(3,100) ' CURRENT DISTRIBUTION FOR CIRCULAR SLIT CYL.' 

NPHI=95 

PHI1=0. 

PHI2=1. 

WRITE(3,110) NPHI 
WRITE(3,120) PHIl 
WRITE(3,120) PHI2 

SMAX=B*(2*PI-ABS(TP2-TP1)-PI/2)+(B-A)+A*(2*PI-ABS(TP2-TPl)-PI/2) 
DO 33 P=0.,.52,.01 
S=P*SMAX/B+PI/2 
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S3=EXP(i*Kl*COS(TP2+S)) 

S4=(0,0) 

DO 22 N=1,NMAX+1 
NMX=FLOAT(N)-l. 

IF(NMX.EQ.O.) THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

S1=(0,0) 

S2=0. 

Sl=EN*CBl*i**(NMX-l.) * DJ(NMX)/(1.+CB2*NMX**2) 
S2=COS(NMX*(S+TP2)) 

DN(N)=S1*S2 

S4=S4+DN(N) 

22 CONTINUE 

JP=-(S4+S3) 

JPHI=CABS(JP) 

WRITE(3,120) JPHI 
33 CONTINUE 

C PLOTTING ROUTINE 

DO 10 R=.57,l.,.01 
S= (R-.57)*SMAX/A 
S5=EXP(i*K2*COS(TPl-S)) 

S4=(0,0) 

DO 5 N=1,NMAX+1 

NMX«FLOAT(N)-l. 

IF(NMX.EQ.O.) THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

S1=(0,0) 

S2=0. 

Sl=EN*CAl*i**(NMX-l.)*DJ1(NMX)/(1.+CA2*NMX**2) 
S2=COS(NMX*(TPl-S)) 

DN(N)=S1*S2 

S4=S4+DN(N) 

5 CONTINUE 

JP=-(S4+S5) 

JPHI=CABS(JP) 

WRITE(3,120) JPHI 
10 CONTINUE 

CLOSE(3) 

100 FORMAT(A) 

110 FORMAT(I5) 

120 FORMAT(E12.3) 

STOP 

END 
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APPENDIX E. COMPUTER PROGRAM FOR /o/X versus 0 FOR SLIT CIRCULAR 
CYLINDERS 


C THIS PROGRAM COMPUTES THE SCATTERING WIDTH SQRT(SIGMA/LAMBDA) 
C FOR A CIRCULAR HOLLOW CYLINDER WITH A SLIT OF ANY SIZE FOR THE 
C TRANSVERSE ELECTRIC POLARIZATION CASE USING THE ON - SURFACE 
C RADIATION CONDITION APPROACH. 

REAL*8 J(0:200),Y(0;200),DY(0:200),DJ(0:200),K1D,K2D 
REAL*8 Jl(0:200),Y1(0:200),DJ1(0:200),DY1(0:200),H1 
REAL EN,NMX,R,U,UU,H,P,PI,PHI,FAC,F(99),P1 
INTEGER M,N,NT,NN,NMAX 

REAL TA,TB,K,K1,K2,TP1,TP2,FIL0N,FR,B,A,L 

COMPLEX P2,P3,i,Sl,Tl,T2,T3,T4,CBl,CB2,CAl,CA2 

COMPLEX C1,C2,C3,C4,C5,C6,C7,C8,CN2,CN3,CN5,CN6,CN7,CN8 

COMPLEX DN(999),S2,S3,S4,S5,CNl,CN4 

EXTERNAL FILON 

OPEN ( 16,FILE='DAT.6') 

WRITE(*,*) 'ENTER THE FREQUENCY(REAL) IN MHZ' 

READ(*,*) FR 

WRITE(*,*) ' ENTER THE FIRST SLIT ANGLE(REAL) IN DEGREES' 

READ(*,*) TA 

WRITE(*,*) 'ENTER SECOND SLIT ANGLE(REAL) IN DEGREES' 
READ(*,*) TB 

WRITE(*,*) 'ENTER NMAX(INTEGER)' 

READ(*,*)NMAX 
TPl=TA/57.29578 
TP2=TB/57.29578 
PI=3.1415927 
DTR=PI/180. 
i=(0,l) 

K=PI*FR/150. 

L=300./FR 

B=.6*L 

A=.5*L 

K1=K*B 

K1D=DBLE(K1) 

K2=K*A 

K2D=DBLE(K2) 

M=98 

H=ABS(TP2-TP1)/FLOAT(M) 

CB1= (8 . *Kl**2-8 . *i*Kl) / (8 . *K1**2-12 . *i*Kl--3 .) 
CAl=-(8.*K2**2-8.*i*K2)/(8.*K2**2-12.*i*K2-3.) 

CB2=-4./(8.*K1**2-12.*i*Kl-3.) 
CA2=-4./(8.*K2**2-12.*i*K2-3.) 

CALL BES(NMAX,K1D,J,Y,DJ,DY) 

CALL BES(NMAX,K2D,J1,Y1,DJ1,DY1) 

C PLOTTING ROUTINE 

WRITE(16,960) 'TE CASE PEC SLIT CYLINDER SCATTERING WIDTH ' 
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NPHI=91 
PHI1=0. 

PHI2=180. 

WRITE(16,970) NPHI 
WRITE(16,980) PHIl 
WRITE(16,980) PHI2 
DO 50 P =0.,180.,2. 

PHI=DTR*P 
IF(P.EQ.O.) THEN 
Tl=i*K*(A-B)*COS(TP2) 

T3=i*K*(A-B)*COS(TPl) 

ELSE 

Tl=COS(TP2)/(COS(TP2)-COS(TP2+PHI))*(EXP(i*K2*(COS(TP2)-COS(TP2+ 
+ PHI)))-EXP(i*Kl*(COS(TP2)-COS(TP2+PHI)))) 

T3=COS(TPl)/(COS(TPl)-COS(TPl+PHI))*(EXP(i*K2*(COS(TPl)-COS(TP1+ 
+ PHI)))-EXP(i*Kl*(COS(TPl)-COS(TPl+PHI)))) 

ENDIF 

T2=EXP((-i)*K2*COS(TP2+PHI))-EXP((-i)*Kl*COS(TP2+PHI)) 

T4=EXP((-i)*K2*COS(TPl+PHI))-EXP((-i)*K1*C0S(TPl+PHI)) 

C EVALUATING 1ST FUNCTION F(U) 

DO 100 U =1,M+1 

UU=TP1+(H*(U-1.)) 

F(U)=COS(Kl*COS(UU+PHI)) 

CONTINUE 
S1=(0,0) 

P1=0. 

DO 200 N=1,NMAX+1 
NMX=FLOAT(N)-l. 

IF(NMX.EQ.O.) THEN 
EN=1. 

ELSE 
EN=2. 

ENDIF 
P2=(0,0) 

PI = FILON(F,NMX,TPl,TP2,M+1,99,1) 

WRITE(*,*) NMX,PI 
P2=EN*i**NMX*DJ(NMX) 

WRITE(*,*) P2 
P3=P1*P2 
WRITE(*,*) P3 
SI = SI + P3 
CONTINUE 
C1=S1 

C EVALUATING 2ND FUNCTION F(U) (CN2=Cl-jC2) & CN4 


100 


200 
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DO 250 U=1,M+1 

UU=TP1+(H*(U-1.)) 

F(U)=SIN(K1*C0S(UU+PHI)) 

250 CONTINUE 

S1=(0,0) 

S2=(0,0) 

P1=0. 

DO 300 N=1,NMAX+1 
NMX=FLOAT(N)-l. 

S3=(0,0) 

S4=(0,0) 

P2=(0,0) 

IF(NMX.EQ.O.) THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

Pl=FILON(F,NMX,TPl,TP2,M+1,99,1) 

S3=(-1.)**NMX*J1(NMX)*DJ1(NMX) 

S4=(-l.)**NMX*i**(-l.)*CA1*DJ1(NMX)**2/(1.+CA2*NMX**2) 

DN(N)=2.*PI*EN*(S3-S4)*COS(NMX*PHI) 
P2=EN*i**NMX*DJ(NMX) 

P3=P1*P2 

S1=S1+P3 

S2=S2+DN(N) 

300 CONTINUE 

CN4=S2 

C2=S1 

CN2*Cl«i*C2 

C EVALUATING 3RD FUNCTION F(U) & CN7 

DO 350 U=1,M+1 

UU=TP1+(H*(U-1.)) 

F(U)=COS(UU+PHI)*COS(Kl*COS(UU+PHI)) 

350 CONTINUE 

S1=(0,0) 

S2=(0,0) 

P1=0. 

DO 400 N=1,NMAX+1 
NMX=FLOAT(N)-l. 

S3=(0,0) 

S4=(0,0) 

S5=(0,0) 

P2=(0,0) 

IF(NMX.EQ.O.) THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

Pl=FILON(F,NMX,TPl,TP2,M+1,99,1) 

P2=EN*CBl*i**(-NMX-1.)*(-1.)**NMX*DJ(NMX)/(1.+CB2*NMX**2) 
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S3=C0S(NMX*TP2)*T2 

S4=(CBl*i**(-NMX-1.)*DJ(NMX))/(1.+CB2 *NMX**2) 

S5=(CAl*i**(-NMX-1.)*DJ1(NMX))/(1.+CA2*NMX**2) 

DN(N)=EN/2.*(“!.)**NMX*(S4+S5)*S3 

P3=P1*P2 

S1=S1+P3 

S2=S2+DN(N) 

400 CONTINUE 

C3=S1 

CN7=S2 

C EVALUATING 4RTH FUNCTION (CN3=C3-jC4) & CN8 

DO 450 U=1,M+1 

UU=TP1+(H*(U-1.)) 

F(U)=COS(UU+PHI)*SIN(Kl*COS(UU+PHI)) 

450 CONTINUE 

S1=(0,0) 

S2=(0,0) 

P1=0. 

DO 500 N=1,NMAX+1 
NMX=FLOAT(N)-l. 

S3=(0,0) 

S4=(0,0) 

S5=(0,0) 

P2=(0,0) 

IF(NMX.EQ.O.) THEN 
EN*1. 

ELSE 

EN=2. 

ENDIF 

Pl=FILON(F,NMX,TP1,TP2,M+1,99,1) 

P2=EN*CBl*i**(-NMX-1.)*(-!.)**NMX*DJ(NMX)/(1.+CB2*NMX**2) 
P3=P1*P2 

S3=COS(NMX*TP1)*T4 

S4=(CBl*i**(-NMX-l.)*DJ(NMX))/(l.+CB2*NMX**2) 
S5=(CAl*i**(-NMX-l.)*DJl(NMX))/(l.+CA2*NMX**2) 

DN(N)=EN/2.*(-1.)**NMX*(S4+S5)*S3 

S1=S1+P3 

S2=S2+DN(N) 

500 CONTINUE 

CN8=S2 
C4=S1 

CN3=C3-i*C4 


C EVALUATING 5TH FUNCTION 

DO 550 U=1,M+1 

UU=TP1+(H*(U-1.)) 
F(U)=COS(K2*COS(UU+PHI)) 
550 CONTINUE 

S1=(0,0) 

P1=0. 
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DO 600 N=1,NMAX+1 
NMX=FLOAT(N)-l. 

IF(NMX.EQ.O.) THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

P2=(0,0) 

Pl=FILON(F,NMX,TPl,TP2,M+1,99,1) 
P2=EN*i**NMX*DJl(NMX) 

P3=P1*P2 
S1=S1+P3 
600 CONTINUE 

C5=S1 

C EVALUATING 6TH FUNCTION (CN5=C5-jC6) 

DO 650 U=1,M+1 

UU=TP1+(H*(U-1.)) 

F(U)=SIN(K2*COS(UU+PHI)) 

650 CONTINUE 

S1=(0,0) 

P1=0. 

DO 700 N=1,NMAX+1 
NMX=FLOAT(N)-l. 

IF(NMX.EQ.O.) THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

P2=(0,0) 

Pl=FILON(F,NMX,TPl,TP2,M+1,99,1) 
P2=EN*i**NMX*DJl(NMX) 

P3=P1*P2 
S1=S1+P3 
700 CONTINUE 

C6=S1 

CN5=C5+i*C6 

C EVALUATING THE 7TH FUNCTION & CNl 
DO 750 U=1,M+1 

UU=TP1+(H*(U-1.)) 

F(U)“COS(UU+PHI)*COS(K2*COS(UU+PHI)) 
750 CONTINUE 

S1=(0,0) 

S2=(0,0) 

P1=0. 

DO 800 N=1,NMAX+1 
NMX=FLOAT(N)-1. 

S3=(0,0) 

S4=(0,0) 

IF(NMX.EQ.O.) THEN 
EN=1. 
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ELSE 

EN=2. 

ENDIF 

P2=(0,0) 

P1=FIL0N(F,NMX,TPl,TP2,M+1,99,1) 

P2=EN*CAl*i**(-NMX-1.)*(“1•)**NMX*DJ1(NMX)/(1.+CA2*NMX**2) 
S3=J(NMX)*DJ(NMX) 

S4=i**(-1)*CB1*DJ(NMX)**2/(1+CB2*NMX**2) 

DN(N)=2.*PI*EN*(S3+S4)*COS(NMX*PHI) 

P3=P1*P2 

S1=S1+P3 

S2=S2+DN(N) 

800 CONTINUE 

C7=S1 
CN1=S2 

C EVALUATING THE 8TH FUNCTION (CN6=C7-jC8) 

DO 850 U=1,M+1 

UU=TP1+(H*(U-1.)) 

F(U)=COS(UU+PHI)*SIN(K2*COS(UU+PHI)) 

850 CONTINUE 

S1=(0,0) 

P1=0. 

DO 900 N=1,NMAX+1 
NMX=FLOAT(N)-l. 

IF(NMX.EQ.O.) THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

P2=(0,0) 

Pl=FILON(F,NMX,TP1,TP2,M+1,99,1) 

P2=EN*CAl*i**(-NMX-1.)*(-1.)**NMX*DJ1(NMX)/(1.+CA2*NMX**2) 
P3=P1*P2 
S1=S1+P3 
900 CONTINUE 

C8=S1 

CN6=C7+i*C8 

SIGMA=1./(8.*PI)*CABS(-K1*(CNl-CN2+i*CN3)-K2*(CN4-CN5+i*CN6) 
+ -T1-T3+CN7+CN8)**2 

SIG=(SIGMA)**.5 
WRITE(16,980) SIG 


50 

CONTINUE 


CLOSE(16) 

960 

FORMAT(A) 

970 

FORMAT(15) 

980 

FORMAT(E12.3) 


STOP 


END 
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APPENDIX F. PROGRAM TO COMPUTE O/wb versus KB FOR SLIT 
CIRCULAR CYLINDERS 

C THIS PROGRAM COMPUTES THE SCATTERING WIDTH/PI*B FOR A CIRCULAR 
C CIRCULAR HOLLOW CYLINDER WITH A SLIT OF ANY SIZE FOR THE 
C TRANSVERSE ELECTRIC POLARIZATION CASE USING THE ON - SURFACE 
C RADIATION CONDITION APPROACH . 

REAL*8 J(0:200),Y(0:200),DY(0:200),DJ(0:200),K1D,K2D 
REAL*8 Jl(0:200),Y1(0:200),DJ1(0:200),DY1(C:200),H1 
REAL EN,NMX,R,U,UU,H,P,PI,PHI,F(99),P1 
INTEGER M,N,NT,NMAX 

REAL TA,TB,K,K1,K2,TP1,TP2,FIL0N,FR,B,A,L 

COMPLEX P2,P3,i,Sl,Tl,T2,T3,T4,CBl,CB2,CAl,CA2 

COMPLEX Cl,C2,C3,C4,C5,C6,C7,C8,CN2,CN3,CN5,CN6,CN7,CN8 

COMPLEX DN(99),S2,S3,S4,S5,CN1,CN4 

EXTERNAL FILON 

OPEN(3,FILE='DAT.71•) 

WRITE(*,*) 'ENTER THE PHI ANGLE(REAL) IN DEGREES' 

READ(*,*) P 

WRITE(*,*) 'ENTER THE OUTER RADIUS(REAL) IN METERS' 
READ(*,*) B 

WRITE(*,*) 'ENTER THE INNER RADIUS(REAL) IN METERS’ 
READ(*,*) A 

WRITE(*,*) ' ENTER THE FIRST SLIT ANGLE(REAL) IN DEGREES' 

READ(*,*) TA 

WRITE(*,*) 'ENTER SECOND SLIT ANGLE(REAL) IN DEGREES' 
READ(*,*) TB 

WRITE(*,*) 'ENTER NMAX(INTEGER)' 

READ(*,*)NMAX 
TPl=TA/57.29578 
TP2=TB/57.29578 
PI=3.1415927 
DTR=PI/180. 

PHI=P*DTR 

i=(0,l) 

M=98 

H=ABS(TP2-TP1)/FLOAT(M) 

C PLOTTING ROUTINE 

WRITE(3,960) 'TE CASE PEC SLIT CYLINDER SCATTERING WIDTH ' 

NPHI=59 

PHI1=0.1 

PHI2=6. 

WRITE(3,970) NPHI 
WRITE(3,980) PHIl 
WRITE(3,980) PHI2 
DO 50 R =0.1,6.,.1 
K1=R*B 
K2=R*A 
K1D=DBLE(K1) 

K2D=DBLE(K2) 

CB1=(8.*Kl**2-8.*i*Kl)/(8.*K1**2-12.*i*Kl-3.) 
CAl=-(8.*K2**2-8.*i*K2)/(8.*K2**2-12.*i*K2-3.) 
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CB2=-4./(8.*K1**2-12.*i*Kl-3.) 

CA2=-4./(8.*K2**2-12.*i*K2-3.) 

CALL BES(NMAX,K1D,J,Y,DJ,DY) 

CALL BES(NMAX,K2D,J1,Y1,DJ1,DY1) 

IF(P.EQ.O.) THEN 
Tl=i*K*(A-B)*COS(TP2) 

T3=i*K*(A-B)*COS(TPl) 

ELSE 

T1=C0S(TP2)/(COS(TP2)-COS(TP2+PHI))*(EXP(i*K2*(COS(TP2)-COS(TP2+ 
+ PHI)))-EXP(i*Kl*(COS(TP2)-COS(TP2+PHI)))) 


T3=COS(TPl)/(COS(TPl)-COS(TPl+PHI))*(EXP(i*K2*(COS(TPl)-COS(TP1+ 
+ PHI)))-EXP(i*Kl*(COS(TPl)-COS(TPl+PHI)))) 

ENDIF 

T2=EXP((-i)*K2*COS(TP2+PHI))-EXP((-i)*Kl*COS(TP2+PHI)) 

T4=EXP((-i)*K2*COS(TPl+PHI))-EXP((-i)*Kl*COS(TPl+PHI)) 

C EVALUATING 1ST FUNCTION F(U) 

DO 100 U =1,M+1 

UU=TP1+(H*(U-1.)) 

F(U)=COS(Kl*COS(UU+PHI)) 

100 CONTINUE 

SI=(0,0) 

P1=0. 

DO 200 N=1,NMAX+1 

NMX=FL0AT(N)-1. 

IF(NMX.EQ.0.) THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

P2=(0,0) 

PI = FILON(F,NMX,TPl,TP2,M+1,99,1) 

WRITE(*,*) NMX,P1 
P2=EN*i**NMX*J(NMX) 

H1=J(NMX) 

WRITE(*,*) P2 
P3=P1*P2 
WRITE(*,*) P3 
SI = SI + P3 
200 CONTINUE 

C1=S1 

C EVALUATING 2ND FUNCTION F(U) (CN2=Cl-jC2) & CN4 
DO 250 U=1,M+1 

UU=TP1+(H*(U-1.)) 

F(U)=SIN(Kl*COS(UU+PKI)) 
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250 CONTINUE 

S1=(0,0) 

S2=(0,0) 

P1=0. 

DO 300 N=1,NMAX+1 

NMX=FLOAT(N)-l. 

S3=(0,0) 

S4=(0,0) 

P2=(0,0) 

IF(NMX.EQ.O.) THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

Pl=FILON(F,NMX,TPl,TP2,M+1,99,1) 

S3=(-1.)**NMX*J1(NMX)*DJ1(NMX) 

S4=(-l.)**NMX*i**(-l.)*CA1*DJ1(NMX)**2/(l.+CA2*NMX**2) 
DN(N)=2.*PI*EN*(S3-S4)*COS(NMX*PHI) 

P2=EN*i**NMX*J(NMX) 

P3=P1*P2 

S1=S1+P3 

S2=S2+DN(N) 

300 CONTINUE 

CN4=S2 

C2=S1 

CN2=Cl-i*C2 

C EVALUATING 3RD FUNCTION F(U) & CN7 
DO 350 U=1,M+1 

UU=TP1+(H*(U-1.)) 

F(U)=COS(UU+PHI)*COS(Kl*COS(UU+PHI)) 

350 CONTINUE 

S1=(0,0) 

S2=(0,0) 

P1=0. 

DO 400 N=1,NMAX+1 
NMX=FLOAT(N)-l. 

S3=(0,0) 

S4=(0,0) 

S5=(0,0) 

P2=(0,0) 

IF(NMX.EQ.0.) THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

Pl=FILON(F,NMX,TP1,TP2,M+1,99,1) 

P2=EN*CBl*i**(-NMX-1.)*(-1•)**NMX*DJ(NMX)/(1.+CB2*NMX**2) 
S3=COS(NMX*TP2)*T2 

S4=(CBl*i**(-NMX-l.)*DJ(NMX))/(l.+CB2*NMX**2) 
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S5=(CAl*l**(-NMX-1.)*DJ1(NMX))/(1.+CA2*NMX**2) 

DN(N)=EN/2.*(-1.)**NMX*(S4+S5)*S3 

P3=P1*P2 

S1=S1+P3 

S2=S2+DN(N) 

400 CONTINUE 

C3=S1 

CN7=S2 

C EVALUATING 4RTH FUNCTION (CN3*C3-jC4) 6 CN8 
DO 450 U=1,M+1 

UU=TP1+(H*(U-1.)) 

F(U)=COS(UU+PHI)*SIN(Kl*COS(UU+PHI)) 

450 CONTINUE 

S1=(0,0) 

S2=(0,0) 

P1=0. 

DO 500 N=1,NMAX+1 
NMX=FLOAT(N)-1. 

S3=(0,0) 

S4=(0,0) 

S5=(0,0) 

P2=(0,0) 

IF(NMX.EQ.0.) THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

P1=FIL0N(F,NMX,TP1,TP2,M+1,99,1) 

P2=EN*CBl*i**(-NMX-1.)*(-!.)**NMX*DJ(NMX)/(1.+CB2*NMX**2) 
P3=P1*P2 

S3=C0S(NMX*TP1)*T4 

S4=(CBl*i**(-NMX-l.)*DJ(NMX))/(l.+CB2*NMX**2) 

S5=(CAl*i**(-NMX-1.)*DJ1(NMX))/(1.+CA2*NMX**2) 

DN(N)=EN/2.*(-1.)**NMX*(S4+S5)*S3 

S1=S1+P3 

S2=S2+DN(N) 

500 CONTINUE 

CN8=S2 
C4=S1 

CN3=C3-i*C4 

C EVALUATING 5TH FUNCTION 

DO 550 U-1,M+1 

UU=TP1+(H*(U-1.)) 

F(U)=COS(K2*COS(UU+PHI)) 

550 CONTINUE 

S1=(0,0) 

P1=0. 

DO 600 N=1,NMAX+1 
NMX=FLOAT(N)-l. 


59 



IF(NMX.EQ.O.) THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

P2=(0,0) 

P1=FIL0N(F,NMX,TPl,TP2,M+1,99,1) 
P2=EN*i**NMX*Jl(NMX) 

P3=P1*P2 
S1=S1+P3 
600 CONTINUE 

C5=S1 

C EVALUATING 6TH FUNCTION (CN5=C5-jC6) 

DO 650 U=1,M+1 

UU=TP1+(H*(U-1.)) 
F(U)=SIN(K2*COS(UU+PHI)) 

650 CONTINUE 

S1=(0,0) 

P1=0. 

DO 700 N=1,NMAX+1 
NMX=FLOAT(N)-l. 

IF(NMX.EQ.O.) THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

P2=(0,0) 

Pl=FILON(F,NMX,TP1,TP2,M+1,99,1) 
P2=EN*i**NMX*Jl(NMX) 

P3=P1*P2 
S1=S1+P3 
700 CONTINUE 

C6 ^S1 

CN5=C5-i*C6 


C EVALUATING THE 7TH FUNCTION & CNl 
DO 750 U=1,M+1 

UU=TP1+(H*(U-1.)) 

F(U)=COS(UU+PHI)*COS(K2*COS(UU+PHI)) 
750 CONTINUE 

S1=(0,0) 

S2=(0,0) 

P1=0. 

DO 800 N=1,NMAX+1 
NMX=FLOAT(N)-l. 

S3=(0,0) 

S4=(0,0) 

IF(NMX.EQ.O.) THEN 
EN=1. 

ELSE 

EN=2. 
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ENDIF 

P2=(0,0) 

P1=FIL0N(F,NMX,TP1,TP2,M+1,99,1) 

P2=EN*CAl*i**(-NMX-1.)*(-1.)**NMX*DJ1(NMX)/(1.+CA2*NMX**2) 
S3=J(NMX)*DJ(NMX) 

S4=i**(-1)*CB1*DJ(NMX)**2/(1+CB2*NMX**2) 

DN(N)=2.*PI*EN*(S3+S4)*COS(NMX*PHI) 

P3=P1*P2 

S1=S1+P3 

S2=S2+DN(N) 

800 CONTINUE 

C7=S1 
CN1=S2 

EVALUATING THE 8TH FUNCTION (CN6=C7-jC8) 

DO 850 U=1,M+1 

UU=TP1+(H*(U-1.)) 

F(U)=COS(UU+PHI)*SIN(K2*COS(UU+PHI)) 

CONTINUE 
S1=(0,0) 

P1=0. 

DO 900 N=1,NMAX+1 
NMX=FIiOAT (N)-1. 

IF(NMX.EQ.O.) THEN 
EN=1. 

ELSE 
EN=2. 

ENDIF 
P2=(0,0) 

Pl=FILON(F,NMX,TP1,TP2,M+1,99,1) 

P2=EN*CAl*i**(-NMX-1.)*(“!•)**NMX*DJ1(NMX)/(1.+CA2*NMX**2) 
P3=P1*P2 
S1=S1+P3 
CONTINUE 
C8=S1 

CN6=C7-i*C8 

SIGMA=1./(4.*K1*PI)*CABS(-K1*(CNl-CN2+i*CN3)-K2*(CN4-CN5+i*CN6) 

+ -T1-T3+CN7+CN8)**2 

WRITE(3,980) SIGMA 


50 

CONTINUE 


CLOSE(3) 

960 

FORMAT(A) 

970 

FORMAT(15) 

980 

FORMAT(E12.3) 


STOP 


END 


C 


850 


900 
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APPENDIX 6. PROGRAM TO COMPUTE THE BESSEL FUNCTIONS 


SUBROUTINE BES(N,X,J,Y,DJ,DY) 

THIS SUBROUTINE WAS OBTAINED DURING EC 4600 FROM PROFESSOR 
MICHAEL A. MORGAN. 

Double precision calculation of ordinary Bessel functions, 
Jn(X) 

and Yn(X), and their first derivative, DJ and DY, for integer 

order "n” from n=0 to N with real argument X. 

REAL*8 J(0;200),Y(0:200),DJ(0:200),DY(0;200),SCALE,JTEMP2,X 
REAL*8 SCLFAC,A,B,C,D,E,F,PI,JTEMP,JTEMPl 
PI=3.14159265359D0 
IF (X.EQ.O.ODOO) THEN 
X = 0.0 BOUNDARY CASE 
IF (N.EQ.l) THEN 

J(l) = O.ODOO 
DJ(1) = 0.5D00 

ELSE 

DO 5, I = N, 2, -1 
J(I) = O.ODOO 
DJ(I) = O.ODOO 
5 CONTINUE 

J(l) = O.ODOO 
DJ(1) = 0.5D00 
ENDIF 

J(0) = l.ODOO 
DJ(0) = O.ODOO 
Y(N) = -l.OD-300 
DY(N) = 1.0D300 
ELSEIF (N.EQ.O) THEN 

POLYNOMIAL EXPANSION ONLY FOR N = 0 
CALL BES0(X,J,Y,PI,DJ,DY) 

ELSE 

RECURSION FOR ALL OTHER CASES 
Y IS A FORWARD RECURSION 
CALL BES0(X,J,Y,PI,DJ,DY) 

Y(l) = -DY(0) 

DY(1)=Y(0) - Y(l)/X 
IF (N.EQ.l) GO TO 20 
DO 10, I = 0, N-2 

Y(H-2) = (2.0D00*(I+1)/X)*Y(I+1) - Y(I) 

DY(I+2) = Y(I+1) - ((I+2)/X)*Y(I+2) 

10 CONTINUE 

J IS A REVERSE RECURSION BASED ON A PAIR OF BESSEL FUNCTION 
POINTS DERIVED FROM A TRUNCATED POWER SERIES EXPANSION. THE 
RECURSION IS THEN SCALED TO A KNOWN VALUE, J1(X). 

20 SCALE = -DJ(0) 

NSAVE = N 
IF (X.LE.N) THEN 
N = 5*N+50 
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0+0 o + 


GOTO 25 

ENDIF 

N = IDNINT(N + X*X + 0.5D00) 

C 

25 A = 1.0DOO/DFLOAT(N+1) 

B = 1.0D00/DFLOAT(N+2) 

C = 1.0D00/DFLOAT(N+3) 

D = 1.0D00/DFLOAT(N+4) 

E = 1.0D00/DFLOAT(N+5) 

F = X/2.0D00 
C 

JTEMP = 1-A*F**2+0.5D00*A*B*F**4-(1.0D00/6.0D00)*A*B*C*F**6+ 
(1. ODOO/24 . ODOO) *A*B*C*D*F**8- (1. ODOO/120. ODOO) *A*^-*C*D*E*F**10 
N = N - 1 

A = 1.0DOO/DFLOAT(N+1) 

B = 1.0D00/DFLOAT(N+2) 

C = 1.0D00/DFLOAT(N+3) 

D = 1.0D00/DFLOAT(N+4) 

E = 1.0D00/DFLOAT(N+5) 

F = X/2.0D00 

JTEMPl = 1-A*F**2+0.5D00*A*B*F**4-(1.0D00/6.0D00)*A*B*C*F**6+ 
(1.ODOO/24.ODOO)*A*B*C*D*F**8-(1.ODOO/120.ODOO)*A*B*C*D*E*F**10 
DO 30, I = N+1,2,-1 

JTEMP2 = 2*((I - 1)/X)*JTEMPl - JTEMP 
IF(DABS(JTEMP2).GE.1.0D250) THEN 
JTEMP2 = JTEMP2*1.0D-250 
JTEMPl = JTEMPl*!.OD-250 

ENDIF 

JTEMP = JTEMPl 
JTEMPl = JTEMP2 
IF ((1-2).LE.NSAVE) THEN 
J(I-2) = JTEMP2 

ENDIF 
30 CONTINUE 
C SCALING 

N = NSAVE 

SCLFAC = SCALE/J(1) 

DO 40, I = 0, N 
J(I) = SCLFAC*J(I) 

IF (I.EQ.O) THEN 
GOTO 40 

ENDIF 

IF (ABS(J(I)/J(I-1)).LT.l.OD-50) THEN 
J(I) * J(I)*1.0D250 

ELSEIF (ABS(J(I)/J(I-1)).GT.1.0D50) THEN 
J(I) = J(I)*1.0D-250 

ENDIF 
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DJ(I) = J(I-l) - (I/X)*J(I) 

40 CONTINUE 
ENDIF 
RETURN 
END 
C 

SUBROUTINE BES0{X,J,Y,PI,DJ,DY) 

C FOR ZERO ORDER BESSEL FUNCTIONS ONLY 

DIMENSION J(0:200), Y(0:200), DJ(0:200), DY(0:200) 

DOUBLE PRECISION J, Y, X, PI, DJ, DY, FO, FI, THETAO 
DOUBLE PRECISION THETAl, A 
IF (X.LE.3.0D00) THEN 
A = X/3.0D00 

J(0) = l.ODOO - 2.2499997D00*(A**2) + 

1.2656208D00*(A**4) - 

+0.3163866D00*(A**6) + 0.0444479D00*(A**8) 

0.0039444D00*(A**10) + 

+0.00021D00*(A**12) 

Y(0) = (2.0D00/PI)*DLOG(X/2.0D00)*J(0) + 0.36746691D00 

+ 

+.60559366D00*(A**2) - 0.74350384D00*(A**4) + 

0.25300117D00*(A**6) 

+ - 0.04261214D00*(A**8) + 0.00427916D00*(A**10) 

0.00024846D00* 

+ (A**12) 

DJ(0)=-X*(.5D00-0.56249985D00*(A**2)+0.21093573D00 
+*(A**4)- 0.03954289D00*(A**6) + 0.00443319D00*(A**8) 

0.00031761 

+D00*{A**10) + 0.00001109D00*(A**12)) 

D Y ( 0 ) 

(“1.ODOO/X)*((2.ODOO/PI)*X*DLOG(X/2DOO)*(-1.ODOO* 

+DJ(0))-0.6366198DOO+0.2212091D00*(A**2)+2.1682709D00*(A**4) 

+1.3164827D00*(A**6) + 0.3123951DO0*(A**8) 

0.0400976D00*(A**10) 

+ + 0.0027873D00*(A**12)) 

ELSE 

A = 3.ODOO/X 

FO = .79788456D00 - 0.00000077D00*A 

0.00552740D00*(A**2) 

+-0.00009512D00*(A**3) + 0.00137237DO0*(A**4) 

-0.00072805D00*(A**5) 

++0.00014476D00*(A**6) 

THETAO = X - 0.78539816D00 - 0.04166397D00*A - 

0.00003954 

+D00*(A**2) + 0.00262573D00*(A**3) - 0.00054125D00*(A**4) - 
+0.00029333D00*(A**5) + 0.00013558D00*(A**6) 

J(0) = FO*DCOS(THETAO)/DSQRT(X) 

Y(0) = FO*DSIN(THETAO)/DSQRT(X) 

FI = 0.79788456D0C + 0.00000156D00*A + 0.01659667D00*A*A 
++0.00017105DOO*(A**3) - 0.00249511D00*(A**4) + 0.00113653D00 
+*(A**5) -0.00020033D00*(A**6) 
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THETAl = X - 2.35619449D00 + .12499612D00*A + 0.00005650 
+D00*(A**2) - 0.00637879D00*(A**3) + 0.00074348D00*(A**4) + 

+ 0.00079824D00*(A**5) - 0.O0029166D0OO*(A**6) 

DJ(0) = -F1*DC0S(THETA1)/DSQRT(X) 

DY(0) * -F1*DSIN(THETAl)/DSQRT(X) 

ENDIF 

RETURN 

END 
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APPENDIX H. PROGRAM TO COMPUTE THE INTEGRATION ROUTINES 


B 

This program computes INTEGRAL F(x) cosTx dx if KEY = 1 and 
B A 

Int F (X) sin Tx dx when KEY=2 using Filon's mwthod. F(x) is 
A 

assumed to be given as a table of functional values in the 
array at H equidistant points from A to B, M odd. MAX is the 
dimension of the array declared in the calling program. 

Ref: P. J. Davis and P. Rabinowitz, Methods of Numerical 
Integration, Academic Press, p. 488 


REAL FUNCTION FILON (F, T, A, B, M, MAX, KEY) 

INTEGER KEY, M, N, I, MAX 

REAL F (MAX), AL, BE, GA, SUM, TH, H, S, C, SI, Cl, S2, C2 

REAL FI, F2, SU, SUl, Al, A, B, T, TH2, TH4 

N = M - 1 

H =ABS(B - A) / FLOAT (N) 

TH = T * H 
S = SIN (TH) 

C = COS (TH) 

TH2 = TH * TH 

IF (ABS (TH) .LE. 0.1) THEN 
TH4 = TH2 * TH2 

AL = 2. / 45. * TH * TH2 * (1. - TH2 / 7. + TH4 / 105.) 

BE =2. / 3. +2. / 15. * TH2 - 4. / 105. * TH4 

GA*4. /3. -2. /15. * TH2 + 1. / 210. * TH4 

ELSE 

AL=1. /TH+S*C/ TH2 - 2. * S * S / (TH2 * TH) 

BE=(2. +2. *C*C -4. *S*C/ TH) / TH2 

GA = 4. * (S / TH - C) / TH2 

END IF 

SUM =0.0 

FI = F (1) 

F2 = F (M) 

51 = SIN (A * T) 

52 = SIN (B * T) 

Cl = COS (A * T) 

C2 = COS (B * T) 

Al = A + H 

GO TO (1, 2), KEY 
1 SU = F2 * S2 - FI * SI 

SUl = -0.5 * (F2 * C2 - FI * Cl) 

DO 3 I = 2, N, 2 

SUM = SUM + F (I) * COS (Al * T) 

Al = Al + H 

SUl = SUl + F (I+l) * COS (Al * T) 

3 Al = Al + H 

GO TO 4 
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If) 


2 IF (T .EQ. 0.) THEN 

FILON = 0. 

RETURN 
END IF 

SU * FI * Cl - F2 * C2 

SUl = -0.5 * (F2 * S2 - FI * SI) 

DO 5 I = 2, N, 2 

SUM = SUM + F (I) * SIN (A1 * T) 

A1 = A1 + H 

SUl = SUl + F (I+l) * SIN (A1 * T) 

A1 = A1 + H 

FILON = H * (AL * SU + BE * SUl + GA * SUM) 
RETURN 
END 
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